#define INF     2e10f
#define DIMX     DIM_X
#define DIMY     DIM_Y
#define ENTITYS ENTITYS_NUM
#define DEPTH 30
#define ABS(x) (x>0.0f ? x:-x)
#ifndef MAX
#define MAX(a,b) (a > b ? a : b)
#endif
#ifndef MIN
#define MIN(a,b) (a > b ? b : a)
#endif
#ifndef M_PI
#define M_PI 3.14159265f
#endif
#define INV_PI 0.3183098861837907f
// //////////////////产生随机数///////////////////////////////////////////////////
/*
Part of MWC64X by David Thomas, dt10@imperial.ac.uk
This is provided under BSD, full license is with the main package.
See http://www.doc.ic.ac.uk/~dt10/research
*/
#define f3(f) f.xyz

//材质
typedef struct Mtl
{
    float4 diff;
	float4 spec;
    float4 emit;
} Mtl;

//顶点
typedef struct Vertex
{
    float4 position;
    float4 normal;
} Vertex;
//三角形面片
typedef struct Triangle
{
    Vertex v1;
	Vertex v2;
	Vertex v3;
    //unsigned int mtlIndex;
    //unsigned int padding[3];
} Triangle;
// 包围盒
typedef struct
{
  float4 pos[2];
} Bounds3;
//bvn树节点
typedef struct LinearBVHNode
{
    // 32 bytes
    Bounds3 bounds;
	float4 info;
    // 4 bytes
    //float offset; // primitives (leaf) or second child (interior) offset,nPrimitives
    // 2 bytes
    //float nPrimitives;  // 0 -> interior node
    // 1 byte
    //unsigned char axis;          // interior node: xyz
    //unsigned char pad[9];        // ensure 48 byte total size

} LinearBVHNode;

// 光线交点
typedef struct
{
    bool hit;
    //Ray ray;
    float t;
    float3 pos;
    //float3 texcoord;
    float3 normal;
	unsigned char mtype;
	unsigned char mindex;
    const __global Triangle* object;
} IntersectData;

typedef struct{ uint x; uint c; } mwc64x_state_t;
typedef struct
{
    float3 origin;
    float3 dir;
    float3 invDir;
    int sign[3];
} Ray;

enum{ MWC64X_A = 4294883355U };
enum{ MWC64X_M = 18446383549859758079UL };

void MWC64X_Step(mwc64x_state_t *s)
{
	uint X=s->x, C=s->c;
	
	uint Xn=MWC64X_A*X+C;
	uint carry=(uint)(Xn<C);				// The (Xn<C) will be zero or one for scalar
	uint Cn=mad_hi(MWC64X_A,X,carry);  
	
	s->x=Xn;
	s->c=Cn;
}
//! Return a 32-bit integer in the range [0..2^32)
uint MWC64X_NextUint(mwc64x_state_t *s)
{
	uint res=s->x ^ s->c;
	MWC64X_Step(s);
	return res;
}


// 计算光线和三角形面皮相交(这个函数是错的！）
bool RayTriangle(const Ray* r, const __global Triangle* triangle, IntersectData* isect,unsigned char mtype)
{
    float3 e1 = triangle->v2.position.xyz - triangle->v1.position.xyz;
    float3 e2 = triangle->v3.position.xyz - triangle->v1.position.xyz;
    // Calculate planes normal vector
    float3 pvec = cross(r->dir, e2);
    float det = dot(e1, pvec);

    // Ray is parallel to plane
    if (ABS(det)<1e-6f)
    {
        return false;
    }
    float inv_det = 1.0f / det;
    float3 tvec = r->origin - triangle->v1.position.xyz;
    float u = dot(tvec, pvec) * inv_det;
    if (u < 0.0f || u > 1.0f)
    {
        return false;
    }

    float3 qvec = cross(tvec, e1);
    float v = dot(r->dir, qvec) * inv_det;
    if (v < 0.0f || u + v > 1.0f)
    {
        return false;
    }

    float t = dot(e2, qvec) * inv_det;

    if (t < isect->t)
    {
        isect->hit = true;
        isect->t = t;
        isect->pos = r->origin + r->dir * t;
        isect->object = triangle;
		isect->mtype=mtype;
        isect->normal = normalize(u * triangle->v2.normal.xyz + v * triangle->v3.normal.xyz + (1.0f - u - v) * triangle->v1.normal.xyz);
        //isect->texcoord = u * triangle->v2.texcoord + v * triangle->v3.texcoord + (1.0f - u - v) * triangle->v1.texcoord;
    }

    return true;
}

// 计算光线和三角形面皮相交
bool intersect_triangle(const Ray* r, const __global Triangle* triangle, IntersectData* isect)
{
	float3 v1=triangle->v1.position.xyz;
	int mtype=(int)triangle->v1.position.w;
	float3 v2=triangle->v2.position.xyz;
	int mindex=(int)triangle->v2.position.w;
	float3 v3=triangle->v3.position.xyz;
    float3 e1 = v2 - v1;
    float3 e2 = v3 - v1;
    // Calculate planes normal vector
	float3 O=r->origin;
	float3 D=r->dir;
    float3 N = cross(e1, e2); //平面法向量
    float det = dot(N, D);
    // Ray is parallel to plane
    if (ABS(det)<1e-6f) //平行
    {
        return false;
    }
	// 计算与平面的距离
    float d = dot(v1 - O, N) / det;
    if (d < 0) //平面在射线后头
        return false;
	if (d > isect->t) return false;
	float3 M=O+d*D; //交点
	// 确定是否在三角形内部
	if (dot(N,cross(e1,M-v1))<0) return false;
	if (dot(N,cross(v3-v2,M-v2))<0) return false;
	if (dot(N,cross(v1-v3,M-v3))<0) return false;
    
	isect->hit = true;
	isect->t = d;
	isect->pos = M;
	isect->object = triangle;
	isect->mtype=mtype; //材料类型
	isect->mindex=mindex; //材料类型
	isect->normal = N;//normalize(u * triangle->v2.normal.xyz + v * triangle->v3.normal.xyz + (1.0f - u - v) * triangle->v1.normal.xyz);
        //isect->texcoord = u * triangle->v2.texcoord + v * triangle->v3.texcoord + (1.0f - u - v) * triangle->v1.texcoord;
    return true;
}

// 计算光线和三角形面皮相交
bool intersect_triangle_uv(const Ray* r, const __global Triangle* triangle, IntersectData* isect)
{
	float3 v1=triangle->v1.position.xyz;
	int mtype=(int)triangle->v1.position.w;
	float3 v2=triangle->v2.position.xyz;
	int mindex=(int)triangle->v2.position.w;
	float3 v3=triangle->v3.position.xyz;
    float3 e1 = v2 - v1;
    float3 e2 = v3 - v1;
    // Calculate planes normal vector
	float3 O=r->origin;
	float3 D=r->dir;
    float3 N = cross(e1, e2); //平面法向量
    float det = dot(N, D);
    // Ray is parallel to plane
    if (ABS(det)<1e-7f) //平行
    {
        return false;
    }
	// 计算与平面的距离
    float d = dot(v1 - O, N) / det;
    if (d < 0) //平面在射线后头
        return false;
	if (d > isect->t) return false;
	float3 M=O+d*D; //交点
	// 确定是否在三角形内部
	float e1e1=dot(v2-v1,v2-v1);
	float e2e2=dot(v3-v1,v3-v1);
	float e1e2=dot(v2-v1,v3-v1);
	float invbase=1.0f/(e1e2*e1e2-e1e1*e2e2);
	float3 p=M-v1;
	float pe1=dot(p,v2-v1);
	float pe2=dot(p,v3-v1);
	
	float u=-(pe1*e2e2-pe2*e1e2)*invbase;
	if (u<0 || u>1) return false;
	float v=(pe1*e1e2-pe2*e1e1)*invbase;
	if (v<0 || v>1) return false;
	if (u+v>1) return false;
	isect->hit = true;
	isect->t = d;
	isect->pos = M;
	isect->object = triangle;
	isect->mtype=mtype;
	isect->mindex=mindex;
	isect->normal = normalize(u * triangle->v2.normal.xyz + v * triangle->v3.normal.xyz + (1.0f - u - v) * triangle->v1.normal.xyz);
        //isect->texcoord = u * triangle->v2.texcoord + v * triangle->v3.texcoord + (1.0f - u - v) * triangle->v1.texcoord;
    return true;
}


//与无限平面
inline float intersect_plane(float3 O, float3 D, float3 P,float3 N)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // plane (P, N), or +inf if there is no intersection.
    // O and P are 3D points, D and N (normal) are normalized vectors.
    float denom = dot(D, N);
    if (ABS(denom) < 1e-6)
        return INF;
    float d = dot(P - O, N) / denom;
    if (d < 0)
        return INF;
    return d;
}

//inline float intersect_sphere(float3 O,float3 D,float3 S,float R) (ray, &face, &isect)
bool intersect_sphere(const Ray* ray, const __global Triangle* triangle, IntersectData* isect)
{
    // Return the distance from O to the intersection of the ray (O, D) with the
    // sphere (S, R), or +inf if there is no intersection.
    // O and S are 3D points, D (direction) is a normalized vector, R is a scalar.
	float3 D=ray->dir;
	float3 O=ray->origin;
	float3 S=triangle->v1.normal.xyz;
	float R=triangle->v1.normal.w;
	float3 SO=S-O;
    float b = dot(D, SO);
    float c = dot(SO, SO) - R*R;
    float disc = b * b - c;
    if (disc > 0 && b>0){
		float rt=sqrt(disc);
        float t=((b-rt)>0 && (b+rt)>0)?(b-rt):(b+rt); //同时处理反射和折射
		if (t<isect->t)
		{
			isect->hit = true;
			isect->t = t;
			isect->pos = O+t*D;
			isect->object = triangle;
			isect->mtype=(int)triangle->v1.position.w;
			isect->mindex=(int)triangle->v2.position.w;
			isect->normal =normalize(O+t*D-S);
			return true;
		}
	}
    return false;
}

bool intersect_obj(const Ray* ray, const __global Triangle* triangle, IntersectData* isect)
{
	bool res;
	if (triangle->v3.position.w>0.5f)
		res=intersect_sphere(ray, triangle, isect);
	else 
		res=intersect_triangle_uv(ray, triangle, isect);
	return res;
}

//暴力算法
inline IntersectData find_intersect_object(Ray* ray ,__global const Triangle* faces,__global const LinearBVHNode* bvh,__global float4* debugarr)
{
    //float t = INF;
    IntersectData isect;
    isect.hit = false;
    //isect.ray = *ray;
    isect.t = INF;
    for (int i=0;i< ENTITYS;i++)
	{
		//debugarr[i]=normal[i];
		// debugarr[i].w=i+1;
		//RayTriangle(ray, &faces[i], &isect,0);
		intersect_obj(ray, &faces[i], &isect);
	}
    return isect;
}

//bvh算法
bool RayBounds(const __global Bounds3* bounds, const Ray* ray, float t)
{
    float t0 = max(0.0f, (bounds->pos[ray->sign[0]].x - ray->origin.x) * ray->invDir.x);
    float t1 = min(t, (bounds->pos[1 - ray->sign[0]].x - ray->origin.x) * ray->invDir.x);

    t0 = max(t0, (bounds->pos[ray->sign[1]].y - ray->origin.y) * ray->invDir.y);
    t1 = min(t1, (bounds->pos[1 - ray->sign[1]].y - ray->origin.y) * ray->invDir.y);

    t0 = max(t0, (bounds->pos[ray->sign[2]].z - ray->origin.z) * ray->invDir.z);
    t1 = min(t1, (bounds->pos[1 - ray->sign[2]].z - ray->origin.z) * ray->invDir.z);

    return (t1 >= t0); //=很有必要，例如厚度为零的包含盒和直线相交

}
inline IntersectData find_intersect_object_bvh(Ray* ray ,__global const Triangle* faces,__global const LinearBVHNode* bvh,__global float4* debugarr)
{
    float t = INF;
    IntersectData isect;
    isect.hit = false;
    //isect.ray = *ray;
    isect.t = INF;
    int toVisitOffset = 0, currentNodeIndex = 0;//等于0代表根节点
    int nodesToVisit[64]; //用栈代替递归
	while(true)
	{
		__global LinearBVHNode* node = &bvh[currentNodeIndex];// 从全局内存读入要测试的当前节点
		if (RayBounds(&node->bounds,ray, isect.t)) //光线和包围盒相交
		{
			
			//如果是叶子节点 info.y= nPrimitives  info.x =offset
			if((int)node->info.y>0) 
			{
				//与该叶子节点中的面片测试相交
				for (int i=0;i<(int)node->info.y;i++){
					//intersect_triangle_uv
					intersect_obj(ray, &faces[(int)node->info.x+i], &isect);
				}
				// 然后右节点出栈
				if (toVisitOffset==0){break;} //栈空了，遍历结束，退出
				toVisitOffset-=1;
				currentNodeIndex=nodesToVisit[toVisitOffset];
			}
			else
			{
				// 考虑划分的方向 大概能减少一半时间
				if (ray->sign[(int)node->info.z]) //光从右边射入，先访问右节点
				{
					nodesToVisit[toVisitOffset]=currentNodeIndex+=1; 
					toVisitOffset+=1;
					currentNodeIndex=(int)node->info.x;;//左节点
				}
				else
				{
					//内部节点，入栈右节点，然后访问左节点
					nodesToVisit[toVisitOffset]=(int)node->info.x; //offset就是右节点在bvh数组中的索引
					toVisitOffset+=1;
					currentNodeIndex+=1;//左节点
				}
			}
		}
		else //光线和包围盒不相交，右节点出栈
		{
			if (toVisitOffset==0){break;} //栈空了，遍历结束，退出
			toVisitOffset-=1;
			currentNodeIndex=nodesToVisit[toVisitOffset];
		}
	}
	return isect;
}

// inline int find_shadow(float3 M,float3 N,float3 toL,int obj_idx,__global const float4* pos,__global const float4* normal,__global float4* debugarr)
// {
	// float ll=length(toL);
    // for(int k=0;k< ENTITYS;k++){
        // if (k != obj_idx){
            // float t = intersect(M + N * 0.0001f, normalize(toL), pos[k],normal[k].xyz,debugarr);
            // if(t<=ll) return 1;
		// }
	// }
    // return 0;
// }

inline float3 get_normal(int obj,float3 M ,__global const float4* pos,__global const float4* normal)
{
    // Find normal.
	float4 p=pos[obj];
    if (p.w==INF)
        return normal[obj].xyz;
    else{
		return normalize(M-p.xyz);
	}
}

inline float3 get_color(int obj,float3 M ,__global const float4* pos,__global const float4* color)
{
    // Find normal.
	float4 p=pos[obj];
    if (p.w==INF)
        return (float3)(((int)(M.x*2+2000)%2)==((int)(M.z*2+2000)%2) ? 0.9f : 0.1f);
    else{
		return color[obj].xyz;
	}
}

//计算反射光线
inline float3 reflect(float3 v, float3 n)
{
	return v - 2 * dot(v, n) * n;
}

// 计算折射光线
inline bool refract2(float3 v, float3 n, float ni_over_nt, float3* refracted)
{
	float3 uv = normalize(v);
	float dt = dot(uv, n);
	// 如果从法线负方向入射（例如球体内部的光线）
	if(dt>0){
		dt=-dt;
		n=-n;
		//ni_over_nt=1.0f/ni_over_nt;
	}
	else
		ni_over_nt=1.0f/ni_over_nt;
	float discriminant = 1.0 - ni_over_nt * ni_over_nt * (1 - dt * dt);
	if (discriminant > 0)
	{
		*refracted = ni_over_nt * (uv - n * dt) - n * sqrt(discriminant);
		return true;
	}
	else
		return false;
}

bool refract(const float3  v, const float3 n, float ni_over_nt, float3* refracted) {
    float3 uv =  normalize(v);
    float dt = dot(uv, n);
    float discriminant = 1.0 - ni_over_nt*ni_over_nt*(1-dt*dt);
    if (discriminant > 0) {
        *refracted = ni_over_nt*(uv - n*dt) - n*sqrt(discriminant);
        return true;
    }
    else
        return false;
}
inline float fresnel(float3 I, float3 N, float ior) 
{ 
    float cosi = dot(I,N); 
    float etai = 1.0f, etat = ior; 
    if (cosi > 0) 
	{ etai=ior;etat=1.0f ;} 
    // Compute sini using Snell's law
    float sint = etai/etat * sqrt(1.0f - cosi * cosi); 
    
    float cost = sqrt(1 - sint * sint); 
	cosi = ABS(cosi); 
	float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost)); 
	float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost)); 
	float kr = (Rs * Rs + Rp * Rp) / 2; 
    
    // As a consequence of the conservation of energy, transmittance is given by:
    // kt = 1 - kr;
	return kr;
} 
// 折射系数
inline float schlick(float cosine, float ref_idx) {
    float r0 = (1-ref_idx) / (1+ref_idx);
    r0 = r0*r0;
    return r0 + (1-r0)*pow((1 - cosine),5);
}

// 产生0~1随机数
inline float RNG(mwc64x_state_t *rng)
{
	return (float)MWC64X_NextUint(rng)*2.3283064365386963e-10f; //2.3283064365386963e-10f
}

inline float3 random_in_unit_sphere(mwc64x_state_t *rng) {
    float3 p;
    do {
        p = 2.0f*(float3)(RNG(rng), RNG(rng), RNG(rng)) - (float3)(1.0f,1.0f,1.0f);
    } while (length(p) >= 1.0f);
    return p;
}

inline void savelight(float3 p1,float3 p2,int * counter,__global float4* debugarr,int i,int j)
{
	if(i==300 && j==137){
		debugarr[*counter*2].xyz=p1;
		debugarr[*counter*2+1].xyz=p2;
		*counter+=1;
		debugarr[0].w=*counter;
	}
}



// 相机
float3 random_in_unit_disk(mwc64x_state_t *seed) {
    float3 p;
    do {
        p = 2.0f*((float3)(RNG(seed),RNG(seed),0.0f)) - (float3)(1.0f,1.0f,0.0f);
    } while (dot(p,p) >= 1.0f);
    return p;
}
Ray InitRay(float3 origin, float3 dir)
{
    dir = normalize(dir);
    Ray r;
    r.origin = origin;
    r.dir = dir;
    r.invDir = 1.0f / dir;
    r.sign[0] = r.invDir.x < 0;
    r.sign[1] = r.invDir.y < 0;
    r.sign[2] = r.invDir.z < 0;
    return r;
}

inline Ray camera_get_ray(float3 lookfrom, float3 lookat, float3 vup, float vfov, float aspect, float aperture, float focus_dist,float s,float t,mwc64x_state_t *seed)
 {
	// vfov is top to bottom in degrees
	float lens_radius = aperture / 2;
	float theta = vfov*M_PI/180;
	float half_height = tan(theta/2);
	float half_width = aspect * half_height;
	float3 origin = lookfrom;
	float3 w = normalize(lookfrom - lookat);
	float3 u = normalize(cross(vup, w));
	float3 v = cross(w, u);
	float3 lower_left_corner = origin  - half_width*focus_dist*u -half_height*focus_dist*v - focus_dist*w;
	float3 horizontal = 2*half_width*focus_dist*u;
	float3 vertical = 2*half_height*focus_dist*v;
	//ray get_ray(float s, float t)
	float3 rd = lens_radius*random_in_unit_disk(seed);
	float3 offset = u * rd.x + v * rd.y;
	origin=origin + offset;
	//ray r={origin + offset, lower_left_corner + s*horizontal + t*vertical - origin - offset};
	Ray r=InitRay(origin + offset, normalize(lower_left_corner + s*horizontal + t*vertical - origin - offset));
	return r;
}



__kernel void ray_tracing( __global const Triangle* faces,__global const LinearBVHNode* bvh,
							__global const Mtl* Mtls,
							__global const float4* lights, write_only image2d_t bitmap,
						   __global float4* debugarr
						   )
{
    // map from threadIdx/BlockIdx to pixel position
    int i = get_global_id(0);
    int j = get_global_id(1);
	// __local Mtl mtls[2];
	
    // int offset = i + j * get_global_size(0);
	// if (offset<2){
		// mtls[offset]=Mtls[offset];
	// }
	// barrier(CLK_LOCAL_MEM_FENCE);
	if (i>=DIMX || j>=DIMY) return;
	// 随机采样 抗锯齿
	mwc64x_state_t seed={i*9000+19797,j*1040+1432};
	RNG(&seed);
    //float3 O =(float3)(ox,oy,oz); //眼睛的位置
	
	int counter=0;//光线计数
	int obj_idx=0;// 与视线相交的实体的编号

	float3 color_ij=(float3)(0,0,0);
	if (i==0 && j==0)
	{
		debugarr[0]=faces[0].v1.position;
		debugarr[1]=faces[0].v1.normal;
		debugarr[2]=faces[0].v2.position;
		debugarr[3]=faces[0].v2.normal;
		debugarr[4]=faces[0].v3.position;
		debugarr[5]=faces[0].v3.normal;
		
		struct LinearBVHNode node =bvh[0];
		debugarr[6]=node.bounds.pos[0];
		debugarr[7]=node.bounds.pos[1];
		//debugarr[7].w=999;
		//debugarr[8].xy=(float2)(999,999);//(bvh[0].offset,bvh[0].nPrimitives);
	}
	int ns=1000 ;//采样次数
	//float3 tmp=(float3)(0.0f);
	for (int s=0;s<ns;s++){
		float reflection=1.0f;
		float3 tmp=(float3)(1.0f);
		// 相机
		float3 lookfrom=(float3)(LOOK_AT_X,LOOK_AT_Y,LOOK_AT_Z);
		float3 lookat=(float3)(LOOK_TO_X,LOOK_TO_Y,LOOK_TO_Z);
		float dist_to_focus = 10;//length(lookfrom-lookat);
		float aperture = 0.0;
		float u = (float)(i + RNG(&seed)) / (float)(DIMX);
        float v = (float)(j + RNG(&seed)) / (float)(DIMY);
		Ray ray=camera_get_ray(lookfrom,lookat, (float3)(NORMAL_UP_X,NORMAL_UP_Y,NORMAL_UP_Z),40,(float)DIMX/(float)DIMY,aperture,dist_to_focus,u,v,&seed);
		float3 D=ray.dir;// 视线的方向
		float3 ORG=ray.origin;
		for (int n=0;n<DEPTH;n++){
			//float dt= find_intersect_object(ORG, D,&obj_idx,pos,normal,debugarr);
			IntersectData insect;
			if (ENTITYS<50)
			{
				insect=find_intersect_object(&ray ,faces,bvh,debugarr);
			}
			else {
				insect=find_intersect_object_bvh(&ray ,faces,bvh,debugarr);
			}
			if( insect.t==INF || n==DEPTH-1)
			{// 不相交直接退出
				float t = 0.5f*(D.y + 1.0f);
				if (n==DEPTH-1)
					tmp*=0;// 光线被困住，全黑
				else
					tmp*=0;//((1.0f-t)*(float3)(1.0f, 1.0f, 1.0f) + t*(float3)(0.5f, 0.7f, 1.0f));//背景色
				break;
			}
			else
			{
			// 交点
			float3 M = insect.pos;
			// 保存视线
			//savelight(ORG,M,&counter,debugarr, i, j);
			// 法线
			float3 N =insect.normal;
			// // 指向光源
			// float3 L=lights[0].xyz;
			// float3 toL =L - M;
			// // 指向眼睛
			// float3 toO = normalize(O - M);
			// // 判断光源是否被挡住
			// int shadowed=find_shadow(M,N,toL,obj_idx,pos,normal,debugarr);
			// toL=normalize(toL);
			// // 材质参数
			//insect.mtype=insect.mtype ;//(N.y>0.99f)?0:1;
			
			Mtl mat=Mtls[insect.mindex];
			//insect.mtype=1;
			// 颜色属性  faceratio着色方法
			//float3 qp=(float3)(((int)(M.x/2.0f+2000)%2)==((int)(M.z/2.f+2000)%2) ? 0.9f : 0.1f);
			//float3 mycolor=insect.mindex==1?(float3)(0.9f,0.9f,0.9f):qp;//(float3)(0.9f,0.6f,0.6f);;////MAX(0.0f,-dot(N,D))*(float3)(0.9f,0.5f,0.3f); //get_color(obj_idx, M ,pos,color);
			//insect.mtype=0;
			
			// // 计算环境光
			// temp+=0.05f;
			// if (!shadowed){
			// // 计算漫反射 Lambert shading (diffuse).
				// temp+=mat.y*MAX(dot(N, toL), 0) * get_color(obj_idx,M ,pos,color);
			// // 计算镜面反射 Blinn-Phong shading (specular).
				// temp+=mat.z*pow(MAX(dot(N, normalize(toL + toO)), 0), 50) * lights[1].xyz;
			// }
			// 按材质计算将要追踪的下一条反射光线
			if(insect.mtype==0)
			{// 0代表光源
				tmp*=mat.emit.xyz;
				break; //遇到光源追踪结束
			}
			if(insect.mtype==1)
			{ // 1代表纯漫反射表面
				D=normalize(N+random_in_unit_sphere(&seed));
				// float3 direction;
				// do{
					// direction=random_in_unit_sphere(&seed);
				// }while(dot(direction,N)<0);
				// D=normalize(direction);
				float scattering_pdf=MAX(0.0f,dot(D,N))*INV_PI;
				//float pdf=0.5*INV_PI;//dot(D,N)*INV_PI;
				float pdf=dot(D,N)*INV_PI;
				tmp*=mat.diff.xyz*scattering_pdf/pdf+mat.emit.xyz; //或者其他方法混合颜色
			}
			else if(insect.mtype==2)
			{ 
			// 2代表反射表面
				tmp*=mat.spec.xyz+mat.emit.xyz; //或者其他方法混合颜色
				float fuzz =mat.diff.w; //模糊因子
				D=normalize(reflect(D, N)+fuzz*random_in_unit_sphere(&seed));
			}
			else if( insect.mtype==3)
			{ // 3代表透明表面
				tmp*=mat.spec.xyz+mat.emit.xyz; //或者其他方法混合颜色
				//float3 outward_normal;
				float ref_idx,reflect_prob;
				ref_idx=mat.spec.w;
				float3 refrac; //折射光线
				if(refract2(D,N,ref_idx,&refrac)) //可以折射
					reflect_prob=fresnel(D,  N,  ref_idx);
				else 
					reflect_prob=1.0f;
				if(RNG(&seed)<reflect_prob)
				{
					 D=normalize(reflect(D, N));
					 //tmp*=reflect_prob; //或者其他方法混合颜色
				}
				else
				{
					 D=refrac;
					 //tmp*=(1-reflect_prob);
				}
				
			}
			
			ORG=M+0.0001f*dot(N,D)*N; //反射光线的起点
			ray=InitRay(ORG,D);
			}
		}
		color_ij+=tmp;
	}
	color_ij/=ns;
	color_ij.x=sqrt(color_ij.x);
	color_ij.y=sqrt(color_ij.y);
	color_ij.z=sqrt(color_ij.z);
	write_imageui(bitmap,(int2)(i,j),(uint4)((int)(MIN(color_ij.x,1) * 255),
	(int)(MIN(color_ij.y,1) * 255),(int)(MIN(color_ij.z,1) * 255),255));
    // bitmap[offset*4 + 0] = (int)(MIN(color_ij.x,1) * 255);
    // bitmap[offset*4 + 1] = (int)(MIN(color_ij.y,1) * 255);
    // bitmap[offset*4 + 2] = (int)(MIN(color_ij.z,1)  * 255);
    // bitmap[offset*4 + 3] = 255;
}